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Abstract 

The evolution of open systems, subject to both Hamiltonian and 
dissipative forces, is studied by writing the nm element of the time (t) 
dependent density matrix in the form 

1 A 

a — 1 

The so called "square root factors", the 7(i)'s, are non-square matrices 
and are averaged over A systems (a) of the ensemble. This square-root 
description is exact. Evolution equations are then postulated for the 
7(f) factors, such as to reduce to the Lindblad- type evolution equations 
for the diagonal terms in the density matrix. For the off-diagonal terms 
they differ from the Lindblad-equations. The "square root factors" ^(t) 
are not unique and the equations for the 7(i)'s depend on the specific 
representation chosen. Two criteria can be suggested for fixing the 
choice of 7(i)'s one is simplicity of the resulting equations and the 
other has to do with the reduction of the difference between the j(t) 
formalism and the Lindblad-equations. 

When the method is tested on cases which have been previously 
treated by other methods, our results agree with them. Examples cho- 
sen are (?) molecular systems, such that are either periodically driven 



1 



near level degeneracies, for which we calculate the decoherence occur- 
ring in multiple Landau-Zener transition, or else when undergoing de- 
scent around conical intersections in the potential surfaces, (ii) formal 
dissipative systems with Lindblad-type operators representing either a 
non-Markovian process or a two-state system coupled to bosons. 

Attractive features of the present factorization method are com- 
plete positivity, the no higher than linear increase of the implementa- 
tion effort with the number of states involved and the introduction of 
randomness only at the start of the process. 

Keywords: Decoherence, Lindblad operators, Landau-Zener crossing, coni- 
cal intersections, Non-Markovian processes 

PACS number (s): 03.65.Bz 

1 Introduction 

The present authors have proposed a variational formulation to study the 
time evolution of the density matrix for situations including both Hamilto- 
nian and dissipative processes [1]. This was based on an expression for the 
time (t)-dependent density matrix, originally devised in [2] and developed 
by [3], which was written formally as 



In this work we apply this "square root" or "factorization" method to the 
investigation of quantum trajectories in a decoherent environment. 

In the above, the "square-root factors" are the non-square matrix j(t) 
whose components are written as 



and its hermitean conjugate. The upper index a designates the system in 
the ensemble and the lower index n the state of the system. In principle, 
both labels run over infinite values, but for bookkeeping purposes we let 
a take A values and n take N values. (Also, in the examples worked out 
in this paper we have taken for A and N finite and small integers.) Thus, 
7 the primary quantity in the formalism is a rectangular N x A matrix. 
As derived in [4] and in other texts following von Neumann's method [5], 
the density matrix is obtained as an ensemble average over all systems. A 



p(t) =1 (t). 1 +(t) 



(1) 



(2) 
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general nm component is 

1 A 

Pnm = J J2 « (3) 
a=l 

[Equation (3) is quoted by numerous quantum mechanics and statistical 
physics textbooks. The relationship of the 7's to the wave-functions in 
an ensemble and other properties of the density matrix are given in the 
appendix. In equation (1) the dot -symbol is a shorthand for the averaging 
(a normalized inner product) over a.] 

A formally identical form is obtained for p from an alternative definition 
of the density matrix of an open system, namely by starting with the total 
density matrix of the system (s) + its interactive environment (e) and then 
taking the trace over the degrees of freedom of the environment [6]. Thus 
the "factorization" of the density matrix involves no approximation. 

To show this, we write a state vector of the combined system as 

I* >s+e= Y.9n\n >s \a > e (4) 

na 

in terms of complete sets of the system and of the environment. From this 
we form the density matrix operator 

p s+e = |tf > s+e < *| s+e = 9n9^\® >e \n > s < m\, < f3\e (5) 

namf) 

Tracing over the environment states gives the reduced density matrix for 
the system alone 

p s = Tr e p s+e = J2< A|/5 s+e |A >= 9n9m\ n >s< m\ s (6) 

A nma 

The nm matrix element of the operator p s is: 

Pnm=<n\p s \m>=^2gX ( 7 ) 

a 

which is of the square root form in equation (3) with 7° = \J~Agn- 

1.1 The non uniqueness of square-root factors 

The square-root factors are not unique. This can easily be seen from equa- 
tion (1) . Suppose I know a matrix j(t) and use it to calculate p(t): 

p(t)= 7 (t)-~f+(t) (8) 
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Obviously I can take any unitary matrix U in the ensemble space and use 
it to calculate another matrix y(£) such that: 



Y(t) = 7 (t)u 



(9) 



Or in index notation: 



in{t) = r n {m 



j3 



a 



(10) 



Now if j'(t) is used to calculate another density matrix p'(t) using equation 
(1) we obtain: 



p'{t) = i{t) • 7 +'(t) = -y(t)U • Uh + (t) = l{t)U ■ U^ + {t) = p(t) (11) 



Hence "/(t) = j(t)U is just another representation of the density matrix 
p(t), the use of either j'(t) or j(t) can not be distinguished experimentally 
and has no physical significance. Furthermore we can choose a U matrix 
evolving in time such that U = U(t). The considerations that dictate which 
is the best representation to use are discussed in the following sections. 

2 Evolution in a dissipative environment 
2.1 A density matrix formulation 

At the present the Liouville-von Neumann-Lindblad (LvNL) equations, that 
are linear in the density matrix and ensure its complete positivity, are fre- 
quently employed for the evolution of the density matrix in the presence 
of dissipative processes. These are written (with the over dot representing 
time derivation and the time dependence temporarily suppressed) as: 



showing the Hamiltonian (H) and dissipative (L) processes [7] . (There may 
be several of the latter, in which case each process is labelled by an index and 
the rate equation contains a sum over the processes. For notational simplic- 
ity we consider a single process and do without an indexed L. A theoretical 
development leading to the above equations is found in [8].) The relation 
of the Lindblad equation to stochastic (Ito or Stratonovich) formulation of 
dissipative processes in quantum mechanics was developed in [9] , further elu- 
cidated in [10] and comparisons between various rate formalism were made 
in [11]. Several extensions of the Monte-Carlo (MC) or unravelling formal- 
ism were made to non-Markovian and other processes (e.g., [12]-[14]). For 
opposing view points we refer to [15]- [22]. 



hp = -i[H, p] + (2LpL + - L + Lp - pL + L) 



(12) 
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2.2 Time development of the 7's 

We now postulate time-evolution equations in the "square root" method. 
The form of the equations is rationalized by the facts that they will correctly 
give the Liouville-von Neumann equations for Hamiltonian processes and 
that they have the form of the LvNL equations for the diagonal elements of 
the density matrix p nn , when there also dissipative terms. The off-diagonal 
matrix elements p nm (n 7^ m) will be discussed after introduction of the 
initial conditions. 

= -iH nr7 ? + L^lTLlAlT ) _1 - L* rn L rsl f (13) 

The star denotes the complex conjugate. The summation convention for 
doubly appearing roman indices is used, but h in a subscript means no 
summation over n (and no summation is implied for multiple Greek indices). 
The inverse (t^*)" 1 is not an element of the inverse matrix of 7; rather, it 
is the inverse of an element of 7. For the conjugate variable one has 

H-W = i^Hrn + m^Lnrl^lTLls ~ lTK S Un (14) 

As discussed in previous texts using the square-root method [2, 3], the rate 
equations for p nn follow from combination of equation (13) and equation 
(14) . In fact, using the left hand side of these expressions and carrying out 
the ensemble averaging, one obtains the time-rate of a diagonal element of 
the density matrix, as follows 



A 

A a 



a 

1 



= {i[p,H}+2Lptf - L^Lp- ptfL}nn (15) 

On the right hand side all products of the 7's an be written in terms of den- 
sity matrix elements and the expression yields precisely the nn- component 
of right hand side of equation (12) . Similarly to solutions of the LvNL equa- 
tions, the trace of p is preserved (Trp(t) = at all times) and the positivity 
of any diagonal matrix elements p nn is ensured, since p nn = j n ■ 7* > 0. 

On the other hand, the dissipative part of the rate expressions for the 
non diagonal components of the density matrix has the form 
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— . / A l Unlr n rm n nr] r !ml In Is ^rs^rm ^rn^rsJs !m 

+ L ns ^^*L* nr (7?* ) ~ 1 7m + Wtf-yT^fr&rSS} 

= {i[p, H] - L^Lp - P L^L} nm + {LB™tf) mrh + (LB* hm L% n 

(16) 

in which the B tensor is defined as 

A a 

This contains terms that include the inverse quantities (7") -1 and cannot 
be expressed in terms of the density matrix. This was already noticed in [2]. 
In the case that m = n we obtain f?™ m = p rs . Thus, while the square root 
method is self consistent, it is not fully equivalent with Lindblad equations. 
(A detailed comparison with different methods is given in section 5.) The 
difference can be clearly formulated by inserting the following expression in 
the curly brackets in equation (16) and subtracting the same from the last 
two terms: 

Kn\ E «* Lsm = {2rfpL} nm (18) 

Then the curly brackets becomes 

{i\p,H]+2l)pL-l)Lp- pl)L} nm (19) 
which is the Lindblad expression, while the difference can be written as 

A a 

with the definition that 

D Ts n ° = (in*)' 1 L nslm L *nr ~ L *sn L rm (21) 

The square-root method leads therefore (in general) to different solutions 
than the Lindblad equation. We repeat that equation (16) is not used to 
obtain the off-diagonal density matrix terms; rather, these are calculated 
directly from the square-root factors. 
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2.2.1 Two components 

To complete this subsection we write down explicitly the equations as they 
should appear for a two component system to be used in the following sec- 
tions (but suppressing the system index a): 

hii = -i(#n7i + #1272) 

2 72* 

- \L 2 i\ 71 - ^21-^2272 + L* 12 — -(£1171 + £1272) (22) 

7i 

^72 = -«(#2i7i + #2272) 

- L* 2 Ln7i - |Li 2 | 2 72 + ^21^7(^2171 + £2272) (23) 

72 

In the case that the Lindblad operator L does not have diagonal elements 
those equations can be further simplified: 

l 1 2 

Hii = -i(Fn7i + #1272) - |^2i| 2 7i + \Li 2 \ 2 -^- (24) 

hj 2 = -i(tf 2l7l + tf 2272 )-|Li 2 | 2 72 + |£2i| 2 ^ (25) 

72 

2.3 The uniqueness of the 7 evolution equation 

According to section 1.1 7 is not unique. 7 can be replaced by an equally 
plausible representation 7' such that a unitary matrix U connects the two: 

im=i'£wm ( 26 ) 

Taking the derivative of equation (26) we obtain: 

Tn(t) = i' n P W$(t) + tf{t)Ug{t) (27) 

This can also be written as: 

im = U$(t)W£(t) +rtWZ(t)(U- l (t)f u ] (28) 

Introducing the expression from equation (26) and equation (28) into equa- 
tion (13) we obtain the result: 

KU$(t)[YP(t) + lXit)U]l{t)(U-\t)t] = U$(t)[-iH nrl f - L* rn L rsl f 

+LnsifirLl r Uf{t){ 1 rUT{t))- 1 ) (29) 
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Multiplying by the inverse matrix of Up and rearranging terms we obtain: 



-h^{t)u»(t){u- l {t)f u (30) 

We see that the third term in the right hand side of equation (30) (corre- 
sponding to the second term in the right hand side of equation (13)) is to 
some degree arbitrary. Furthermore, if the matrix U(t) is chosen to be time 
dependent a fourth term can be added without effecting the results. How 
should one choose a matrix U{t)l Some recommendations can be given: 

1. Choose a matrix U(t) to simplify equation (30). 

2. Choose a matrix U{t) to avoid singularity condition that may occur 
in the third term in the right hand side of equation (30). More about 
the problem of singular initial condition in the next subsection. 

3. Choose a matrix U (t) in order to reduce the difference between the 
Lindblad formalism and the factorization matrix formalism. This can 
be done in terms of the difference tensor D™™ a defined in equation 
(21). Inserting the expression from equation (26) into equation (21) 
we obtain the result: 

U rs [ U f3\ — \lh ) L nsl m ^hr ~ L sn L rm 

= (irKHtT'LnsiiTKHmir - KnLr m (31) 

Hence the U(t) can be chosen in order to obtain smaller D™ na tensors. 

For the examples which were worked out in this paper it was found 
that choosing the matrix U(t) as the identity matrix produced both simple 
equations and also good agreement with the Lindblad theory. However, for 
more involved cases a different choice of the matrix U(t) may be needed. 

2.4 Initial conditions 

The presence of the inverse in the dissipative part of the evolution equations 
causes the initial conditions (IC) to be of great importance. This is evident, 
because the rate expressions are singular whenever a component probability 
is zero. We can simplify the treatment by expressing the density matrix at 
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t = in a diagonal form. This can always be done by a suitable transfor- 
mation of the basic states. The cases of pure and mixed systems are then 
easily distinguished. The physical IC are, quite generally, that 



K(0)| = V^ forn = l,...,M 



(32) 



and zero for the rest of the states. Here the initial probabilities p\ 



•n 



M 



(33) 



n=l 



by normalization of the density matrix. In a pure state 



M = 1 



(34) 



and in a mixed state 



M > 1 



(35) 



We have already noted that a zero value of a 7-variable at t = (or at 
any later time), causes a singularity on the right hand side of equation (13) 
and equation (14) . This is overcome by starting the integration at a time 
arbitrarily close to and slightly above t = 0, subject to the IC's given by 



with the phases (pa,n taken to be real, but with no other restriction on them. 
When / 0, the second term in the above equation can be ignored. But if 
the first term is zero, the following term is essential. The correction to this 
term can be shown to be, for small t, of the order \L\ 2 t. It can be checked 
that the above choice of IC ensures that the singularity on the right hand 
side of the evolution equation for any component 7 is exactly cancelled by 
the singular time derivative on the left hand side. (For a historical remark, 
singularities of the solutions in the density matrix at t = were noted early 
on, in section 4 of [24]. The fast "slippage" or initially irregular behavior of 
solutions was recognized in [25].) 

We now assume that all solutions of the rate equations for the 7's cor- 
respond to a physical system in the ensemble. The ensemble averaging has 
therefore to be carried out for all possible choices of the initial phases. 




9 



2.4.1 Short time behavior 



Let us now calculate the time development of the density matrix at short 
times, for an initially pure state, in which 

Pll (t = 0) = 1, pi„(0) = 0, Pnm {0) = , for n + 1 + m (37) 

In this case 

\fp\ = 1 (38) 

and 

= n ^ 1 (39) 
For the case n / 1 equation (36) takes the form 



l_im 7 «(t) = e^pL^t = e^\L n l\V2t (40) 
The density matrix has the following short time behavior: 

limp nm = lim 1 ^ T^W-CW « ^(*)^ E e*^--*--) (41) 
where 

fem(t) = |L m i|\/2t /or ra = l,m > 1 

= 2|L„i||L m i|t forn^l^m (42) 

The last case includes diagonal matrix elements for the initially absent 
states, n = m > 1, but here the phase factor averaging over a (the sys- 
tem labels) in equation (41) gives unity. These matrix elements are thus 
seen to give rise to a non-zero value within a time of the order of t ~ L~ 2 . 
Due to the dissipative mechanism (represented by L) , the system will there- 
fore become mixed beyond this time, so that Tr p 2 < 1. This pure-to-mixed 
transition can also be obtained from calculation of the time derivatives of 
the square root factors. (With a purely Hamiltonian interaction (H), an 
initially pure state for this Hamiltonian will persist to be a pure state at 
later times.) 



3 Decoherence in Molecular Systems 
3.1 Driven Level Crossing 

The changes in the density matrix involving two states during fast level 
crossing by a molecular systems were considered in [26] including also dis- 
sipative forces. It was later remarked in [27] that the "square root operator 
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method of [2]", represents an alternative way of showing how a dissipative 
term in the Hamiltonian can cause decoherence. We now formulate the 
evolution equations in this square-root scheme and solve the resulting equa- 
tions. The solutions for the diagonal term in the density matrix, shown 
graphically below in Figure 1 and similar to those shown in [1], have all 
the features appearing in the earlier treatments [26, 27] based on entirely 
different schemes of solution. 

The level crossing (LC) system is characterized by a Hamiltonian (Hlc) 
and a non - Hermitian dissipative (Alc) part written for the two states in 
the matrix forms, as follows: 

„ / ^Gcos(u;t) J \ . . 

\ J — ^Gcos(o;t) j 

Llc = ^v[\ J) (44) 

G denotes the strength of the periodic driving field; the coefficients J and 
r of the tunnelling and of the dissipative mechanisms are denoted by same 
symbols as in [26] and [27]. 

For a Markovian process the rate equations are now written for the two 
elements (7° , r y%) in the density matrix as 

ihtf = iccosMhf + J72 -ir[7f-|72l 2 /7i*] 

ihi% = -lGcosM)7 2 a + J7i a -ir[ 72 a -|7i a | 2 /7 2 a 1 (45) 

We have already called attention to the divisors 7* on the extreme right, 
characteristic of the factorization formalism for dissipative processes, equa- 
tion (14) . The trace of the density p(= 7 • 7 + ) stays constant during the 
motion, by construction. 

We next solve two equations for the j's, subject to the pure state initial 
conditions (7" | 2 = 1, 72 ~ at t = 0. Then from the solutions we form the 
diagonal matrix elements of the density matrix and finally show the results 
in figure 1. For a beginning, the three lower drawings in the figure arose 
from calculations that were carried out for a density matrix referring to an 
"ensemble" consisting of one system. This means that one has A = 1 and the 
summation over a in equation (3) is trivial. The quantity changed between 
the upper three drawings is the strength T of the dissipative term. As this 
increases, a transition takes place from the slow to the fast decoherence 
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regime. We have already noted the remarkable similarity of the results 
obtained here by the factorization method to those in Figure 1 in [26] and 
in [27], except that for strong dissipation their drawings show only tiny 
oscillations, unlike our third drawing from below. In this, drawn for — = 20, 
after a very steep initial downward slope (not visible in the figure), both 
diagonal density matrix elements oscillate about the asymptotic value of \. 

In this trivial "ensemble" of A = 1 the off-diagonal matrix element is 
simply given by |pi 2 | = 1 71 72 1 = 1 7i 1 1 72 1 = VPni. 1 ~ Pn) 1 - This differs 
from the same quantity calculated and shown in Fig. 2 of [27]. Their curve 
initially oscillates around 0, and at longer times settles down to this value, 
whereas our result tends asymptotically to \. We do not regard this as a 
serious discrepancy, since in a realistic model, in which A is large, the phase 
decoherence will make pu vanish, as argued in section 2.3. 

We have also calculated the density matrices when there is a non-trivial 

. \ 12L 2i7r Zi-K 4i-K 

summation, namely when initial conditions are 7° (0) = e 2 ,e 2 , e 2 ,e 2 
for a = 1,2,3,4 respectively and 72(0) ~ 0, (so that A = 4), instead of 
having only 7" (0) = 1 (and A = 1), as before. Physically, this represents a 
certain type of averaging over the environment. (In more complex systems, 
one would require an averaging in the density matrix over a much larger 
number of states, such that A — > 00.) The resultant density tends now 
to an almost perfectly straight line. This is similar to the graphs shown 
in both [26] and in [27] for strong dissipation and elucidates the practical 
consequence of system-averaging in the density matrix. We have also worked 
out the A = 4 case for the lower two graphs in Figure 1. For these graphs, 
there was hardly any perceptible change from those shown. 

For the physical meaning of these results in the context of molecular level 
crossing, we point to [26] and references therein, while for a more general 
application to decoherence we refer to [27]. 

3.2 Descent Across a Conical Intersection 

Whereas in the previous case the decoherence mechanism was activated by a 
non-Hermitean dissipative force, whose strength was designated by T, there 
are situations where dissipation is intrinsic in the dynamics. A situation of 
this type takes place in the excited state dynamics of a molecular system 
moving on a potential energy surface as it approaches a conical intersection 
(ci). (This is a frequently encountered natural process, and it has been 

Vn + P22 = 1 => I72! = ^/p22 = yi - pTi 
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Figure 1: Evolution of a diagonal element of the density matrix p\\. The 
strength of dissipation increases upwards in the drawings. For all drawings 
the parameters in equation (45) are chosen as G = 25, J = 3, to = 1. Then, 
in the bottom drawing T (the dissipative parameter)= 0, A (the ensemble 
size) =1; in the second drawing (from below): T = 0.05, A = 1; in the third 
drawing: r = 20, A = 1; in top drawing T = 20 (as in previous, but) A = 4. 
The initial downward slope in the top two drawings is too steep to be visible 
and so are tiny fluctuations in the horizontal part of the top drawing. Note 
the smoothing effect of the ensemble averaging, evident from a comparison 
between the top two drawings. 
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claimed that it is basic to many naturally occurring life processes in which 
an electronic state is changed, e.g. photosynthesis [28]. ci have been previ- 
ously studied by numerous researchers, [29] and other papers in that volume 
[30].) We now briefly give the underlying formal background, with a view 
of applying to it our formalism. A schematic illustration of a ci is shown in 
Figure 2. 

We first recall the basic result of von Neumann and Wigner [31] that the 
crossings (points of degeneracy) of potential energy surfaces (=electronic 
energies as functions of the nuclear coordinates) for a polyatomic molecule 
can be generally described in terms of just two nuclear displacement pa- 
rameters. (Following the notation of [32], we shall here denote these by Q\ 
and Q c , the former being a "tuning" and the latter the "coupling" mode 
coordinate. They are actually two linearly independent combinations of 
the nuclear coordinates.) The two surfaces (belonging to the two locally 
adjacent "adiabatic" electronic states) separate from each other near the 
intersection point in a manner that is linear in both coordinates: hence the 
name "conical intersection" (ci). They differ from energy-surface intersec- 
tions that happen during molecular or atomic collisions, in that these may 
(approximately) be treated in a single-dimensional coordinate space, usu- 
ally the separation between the colliding-reacting species. The probability 
of change of the electronic state in such collisions is given by the Landau- 
Zener formula [33], [34], which is at the basis of the subject in the previous 
section. An analogous semi-classical formula for the passage across a ci was 
obtained by Nikitin [35]. Later developments were summarized in [36]. In 
a simplified form the expression of [35] for the asymptotic probability of 
transfer between the two diabatic states in a 2 -dimensional space (x, y) can 
be written as 



K(vl + 



P 2 11 = exp| - * ' fl 1 (46) 

In this formula v x and v y are components of the starting velocities on the 
upper ("2") adiabatic surface, I is the closest distance in the passage to the 
intersection point and K is the strength of the electron-nucleus dynamic 
coupling. 

The time dependent Schrodinger equation for the dynamics of a ci was 
subsequently solved numerically in several papers, especially in [37] - [38], 
which contain references to earlier works. These show that the initial wave 
packet (which is excited in an upper electronic state) has a definite, non-zero 
asymptotic probability to end up in the other electronic state as the wave 
packet traverses the ci. The dynamics bears therefore the irreversibility 
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Figure 2: Schematic drawings of adiabatic potential surfaces exhibiting con- 
ical intersection (ci) The upper and middle drawings show the upper and 
lower sheets as function of the two displacement coordinates Q\ and Q c intro- 
duced in the text. The lower drawing superimposes these two and shows the 
locally conical nature of the intersection at Q\ = —4 (expressed in arbitrary 
units, typically about -0 .1 nm) and Q c = 0. The system starts its trajectory 
on the upper surface as a wave packet centered at (say) (Qi = —10, Q c = 2) 
(shown by a blob, where it belongs almost entirely to the upper electronic 
state) and descends towards and beyond the ci, oscillating to and forth and 
(partly) losing its upper electronic state character. 
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hallmark of a dissipative mechanism, though such mechanism was nowhere 
introduced in this model (unlike that in the previous section). The irre- 
versibility was explained in [38] as due to the essential anharmonicity of the 
dynamics in the Q±, (^-coordinates. The anharmonicity shows up in the 
cusps on the energy surfaces at the degeneracy point. 



3.2.1 A simplified formalism for ci 

It is customary to represent the two bare or "diabatic" electronic states 
(those that are independent of the nuclear coordinates) by the column vec- 
tors 

("1","2") = (( I )'(?)) < 47 ) 

([32]- [39]). Temporarily simplifying the situation of [32] to the case where 
there is only one tuning mode Qi, as well as the coupling mode Q c , we have 
for the two coordinates the following harmonic oscillator Hamiltonian: 



hui , d 2 ~, Tiu c , d 2 



where %u\ and huj c are the quanta of vibrational energies. In the doublet 
representation the nuclear Hamiltonian is written as a scalar or, equivalently, 
as H nuc times the 2x2 unit matrix /. 

The nuclear-electronic interaction part can then be written without loss 
of generality as 

_ / SE + K 1 Q 1 K C Q C \ 
eZ " nuc ~ ^ K C Q C -5E - i^Qi J { y> 

where 25E is the vertical energy offset between the two electronic states and 
K\,K C are the coupling strengths for the two coordinates. 

The Schrodinger equation to be solved for the time dependent vector 
V(Qi,Q c ,t) is 

. h d*(QuQ c ,t) = + Hel _ nuc) mQ 1:Qct) (50) 

subject to the initial condition at t = 0, that ^f(Qi, Q c , 0) is a wave packet 
(say, of gaussian forms in the two coordinates Qi,Q c ) on the upper po- 
tential surface, centered at some position well above the intersection point 
{Qi = — 7^7 j Qc = 0). This type of wave packet can be conveniently formed 



16 



by short duration (compared to the inverse frequencies and lo^ 1 ) opti- 
cal excitation from some lower lying electronic state. The wave packet will 
then descend towards the ci point and beyond. Throughout its passage the 
initially excited (diabatic) electronic state (say, "1") will lose amplitude in 
favor of the other state "2". This loss will be especially intense near the de- 
generacy point. Ultimately, for "long" times of the order of picoseconds, the 
probability will tend to some asymptotic limit (say, Pi j0 o) between and 1. 
(Figures 1-4, in [32]. Actually, in any individual run it will slightly oscillate 
about the asymptotic probability.) This value depends on the parameters of 
the system and also, presumably to a lesser extent, on the initial conditions. 

3.2.2 A Lindblad-type formalism for ci 

In this section we present a reformulation of the previous two state-two mode 
system, such that the irreversible nature of the dynamics is built into the 
model (rather than emerges from the solution). We do this by constructing 
a phenomenological dissipative term in the square-root formalism, based 
on the Hermitean formalism of the previous subsection. This procedure 
will then lead to Lindblad-type terms in the master equations. (Let it be 
emphasized that we are not deriving the Lindblad term from a microscopic 
process, but are formulating the microscopic equations in an irreversible 
setting for the partial, electronic degrees of freedom, having eliminated the 
nuclear coordinates.) 

In the motion of the wave packet on an adiabatic surface the electronic 
state amplitudes will depend on the position of the nuclear coordinates. In 
the square root formalism, and relying on Appendix A.l, we write this as 

l? = lf{QiM (* = 1,2) (51) 

We invert these relations to make the coordinates some functions of the 
7's and then linearize the functional relation to be composed of a classical 
(non-dissipative, real) part and a "dissipative" part which will give rise to a 
Lindblad type term. Specifically, 

Q r (t) = Q? ass (t) + QL{t) (r = l,c) (52) 

for either coordinate and then 

Qf ass (t) = Q° 1 cos(u 1 t) + — sin(wit) (53) 

Qf ass (t) = ^sin(^) (54) 
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These relations are appropriate for a classical motion (for which the initial 
values of the coordinates are (Q?,0) and of the velocities are (vi,v c ) ); they 
should provide a fair enough description of the center of a gaussian wave 
packet. Other choices for the initial values lead to similar results for the 
asymptotic behavior. (But take note that the level crossing probability is 
zero when the initial velocities are zero or when the classical pathway goes 
across the ci. This is evident from Nikitin's expression for a ci level crossing 
probability [35, 36], shown above in equation (46) .) 

The expressions for the Lindblad component of the coordinates are more 
complicated and their rationale will be apparent only later. They are 

n L (f) g [(1 ~ Mlg - A,oo|7 2 T](7 2 a 7r + 7f 7?*) ^ 

VlU K ± |7fl 2 7?7r + l7 2 T7 2 a 7r 1 ' 

nL(f) «T [(1 ~ Pi^Ml - p i.°°l7 2 a l 2 ](7f 7f* ~ 7 2 Q 7 2 a *) f , fil 

VcU K c l7?| 2 7f7r + l7 2 a ! 2 7 2 a 7r 1 ' 

Here Pi )00 and P 2j00 are asymptotic weights or probabilities for the two 
(adiabatic) components. The time dependence of the 7-factors is suppressed 
in these formula, for brevity. Because of the irreversibility of the process, 
one can assume the functions Qi, Q c to be complex. 

Substitution of the nuclear coordinates, as given in equation (52) - 
equation (56) , into equation (49) and the use of this 7-dependent Hamil- 
tonian in the first ("Hamiltonian") term in equation (13) -equation (14) , 
yield after some simplification the following Lindblad form of the rate of 
equations for the 7's 

ihj? = [5E + Kx (Q? cos(wi*) + — sin(^it))] 7 f + K c — sm(uj c t)^ 

_ ^ ((l-Pl,oo)|7f| 2 -Pl,oo|7 2 a | 2 ) (g7) 



7f* 



ihj% = [-5E - K X (Q\ cos(wit) + — sin^t))]^ + K c — sin{u c t)^ 

+ , r ((l-Pl,oo)l7l| 2 -Pl,oo|72| 2 ) (5g) 

7 2 

The Lindblad operator for this set of equations is: 
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Figure 3: Time development of the weight of the diabatic-state, \^/f (t)\ 2 , 
following excitation into state "1" and during a passage across a coni- 
cal intersection with a lower lying state ("2"). Lower curve: a single 
run. Upper curve (displaced upward by 1): ensemble averaging over four 
runs (A = 4, explained in text). The parameters in equation (57) are 
( with frequencies and energies /% in inverse femtosecond units, lengths 
in units of zero point motion amplitudes, typically .05 nanometers ): En- 
ergy offset, 5E = 0.2, Coupling constants: K\ = 0.06, K c = 0.05. Mode 
frequencies: to\ = uj c = 0.4. Initial velocities: v\ = 0.24, v c = 0.32. Dissipa- 
tive strength, V = —0.1; Asymptotic weight: P\ jOC = 0.25 

When at t = the "1" -state is excited optically, the initial conditions are 
7i(0) a = l,7f (0) = f° r an a i one obtains asymptotically, after some short- 
time oscillations and rise of the occupation probability of the non-excited 
"2" -state, 

Pll(oo) = |7l(°o)| 2 = A,oo, P22(oo) = |72(oo)| 2 = ^2,oo = 1 ~ P\,oo (60) 

The behavior of pn as function of time is shown in figure 3. The parameters 
K\, K c , T, u\, uj c are all functions of the parameters in the molecular Hamil- 
tonian (equation (48) , equation (49) ). Qi,v\,v c are defined by the mode 
of excitation, Pi j0 o depends on the molecular Hamiltonian and, as shown in 
section VI of [38], approximately by the potential surfaces. 

It is proper to characterize the present computed case as belonging to the 
moderately strong dissipative case, since the oscillations about the asymp- 
totic value begin after a time of the order of a vibrational period . A closer 
look at the curves shows that the computed mean asymptotic line lies above 
the "input" asymptotic weight Pi j00 = .25 by about 0.03. This is an inter- 
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esting effect, likely to be due to oscillations in the asymptotic probability 
which arise from the Hamiltonian part of the rate equation (57) . (We call 
attention to the discussion in section VI of [32] , concerning the discrepancy 
of about the same magnitude and sense between their classical and com- 
puted asymptotic probabilities.) Trivially, the other diagonal matrix (P22) 
tends to a value that is lower than its asymptotic value, since the two di- 
agonal matrix elements add up to unity. The fluctuations have a roughly 
uniform period of 2it/uj\. 

The lower curve is for a single run (A = 1), the upper curve (vertically 
displaced for clarity by unity) is after "ensemble averaging" over four runs 
(A = 4), which differ from each other by having different phases (0f = 0, 
7r/2, 7r, 37r/2, for a = 1,..,4) in the initial 7f(0). As seen in Figure 3, 
no qualitative difference is observed between the curves for a single system 
and those for the ensemble averaged density matrix, showing that the single 
system fluctuations due to the periodic Hamiltonian term are not averaged 
out by the dissipative term. 

This last finding changes radically when we apply our formalism to a 
more complex case, that was proposed in [32] as representative of a C2H4 
molecular system. There are now two tuning modes (designated by the 
subscripts 1 and 2), as well as different frequencies for the three modes. The 
results for the diagonal density matrix are shown in Figure 4. 

The fluctuations are here considerably more congested than in the pre- 
ceding case, where there was only a single mode-frequency. One sees in 
Figure 3 that the dissipative mechanism does not completely eliminate the 
fluctuations. However, additional computations (not shown here) indicate 
that increasing the dissipation strength T does further slightly smoothen the 
intrinsic fluctuations. On the other hand, in the ensemble averaged probabil- 
ities (the upper curve in Figure 4) the fluctuations due to the quasi-periodic 
Hamiltonian are almost completely washed out. We finally note that the 
computed mean asymptotic value still lies about 0.03 above the nominal 
asymptotic weight. 

4 Some Models of Dissipative Processes 

4.1 Non-Markovian processes with memory 

Following a work by Diosi et al. [40], we consider the situation in which 
the Lindblad operator L is dependent on time and hence has the faculty of 
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Figure 4: Computations in a model for C2H4 (for a single run, the lower 
curve, and with an "A = 4 ensemble averaging" as in the previous figure, 
the upper curve displaced by unity for clarity). The parameters of Table 
I (A) and Figure 1 (a) of [32] are used with added values for initial ve- 
locities and asymptotic probabilities (chosen somewhat arbitrarily), as fol- 
lows: 6E = 0.95, Coupling constants:^ = 0A9,K 2 = -0.27, K c = 0.5, 
mode frequencies: uj\ = 0.36, uji = 0.21, u c = 0.11, velocities: v\ = .6, v 2 = 
0, v c = .8, T = —.1 (in units as in the previous figure); asymptotic weight: 
Pi,oo = 0.25 
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Figure 5: Expectation values of Pauli matrices for a two - level system 
undergoing non -Markovian quantum diffusion (NMQD) for parameters as 
in Fig. 2(b) of [40] (w = 1; L = 1; G = 1; Wl = 1; Gl = G; a = 0). Full line 
: =< a z >; short, broken lines : =< a x >; long - broken lines : =< a y > 



memory. We consider the simple case of a two-level system, in which: 

H = -wa z 

L = y/f{t)a- (61) 
The memory function is given by: 



G, \Ui-2Gil' l I G-\ 

/(*) = " }LJ -^ tanh[-t JG? - 2G ± P + arctanh[ 1 ]] 

2 2 2 sjGl-IGxX 1 

(62) 

The above Hamiltonian and Lindblad operator lead to the equations: 

7? = -^7f-/(i)7? 

7 2 a = W + (63) 

z 72 

The time evolution of the expectation value of the Pauli matrices is plotted 
in Figure 5 in the case that the above differential equations are solved for 
the pure state initial conditions 7f(0) = 7^^,72(0) = independent of 
a. The obtained results are very similar to those of [40]. 
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4.2 Probabilities for a two-state system coupled to a bosonic 
reservoir 



This problem was recently studied, alongside with other instances also in- 
volving memory, in [41, 42]. The case that we numerically solve using the 
square-root formalism possesses an exact solution [shown in Eq. 57 of [41]], 
with which our computed values agree perfectly (Figure 7). 

In the model a degenerate two-state system, |1 > and |2 > (whose Hamil- 
tonian is taken as zero), is coupled to a reservoir of bosons through a time- 
dependent Hamiltonian interaction term. This is given by 

H in t = B(t)\2 >< 1\ + B*(t)\l ><2| (64) 

where B*(t) is a reservoir excitation operator [Eq.(52) in [41]]. The essential 
quantity in the dynamics is the spectral density given by 

1 IG 2 
2ir (u - ojq) 1 + G\ ) 

where I is the interaction strength, G\ is the spectral width and loq the 
characteristic frequency in the reservoir. The analytical solution for the 
probability of the state occupation |2 > that is excited at t = 0, is 

p 22 (t) = e-^[cosh^ + ^sinh^] (66) 



where ft = ^G\ - 2l 2 G 1 . 

In our formalism, we solve numerically the master equations for 7" (i), 
72 {t) and their complex conjugates, employing a Lindblad operator L = 
f(t)\2 >< 1|, with the memory function f(t) given in equation (62) . Initial 
conditions are the same for all 2 a. The result is shown in Figure 6, which 
completely overlaps the exact result from equation (66) . The Monte Carlo 
method of [41, 42] also agrees with the analytic results, but uses 10 7 runs. 
It is not clear to us what is the minimum runs that is required to obtain 
agreement, but the expediency of the square root method is evident. 

In a further figure we also show the expectation value of the three 
angular momentum matrices taken with respect to the state amplitudes 
(li (t)i I2 (*))■ (Figure 7). They are similar to the drawings in Figure 6, but 
with more persistent oscillations. 

2 This means that the ensemble contains only a single system 
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Figure 6: Probability of excited state occupancy P22(t) = 1 72 (*)| 2 against 
normalized time (u>t). The parameters are I = 1, G\ = 0.2 . Results are 
computed by the square root method. They overlap completely the analytic 
expression. 



a.m. 




wt 



Figure 7: Averages of Pauli, angular momentum matrices for the solutions 
against reduced time. < a z >: full lines; < a x >: long broken lines; < a x >: 
dashed lines. Parameters are as in Fig. 6 
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5 Comparison with Other Methods. 



In this section we briefly consider similarities and differences between the 
"square-root method" and some other methods, already mentioned in the 
Introductory section. 

(I) Number of independent rate equations. 

At any instance of time the density matrix, whether by its definition in 
equation (3) or by environment tracing , is composed of N x A independent 
(complex) numbers (N being the nominal number of states and A either the 
number of systems in the ensemble or the total number of states in the 
environment). Thus, in principle, one needs just N x A rate equations to 
obtain the density. The square root method postulates N x A equations. 
However, the roles of N and A are different and we find that frequently one 
can get good approximations (to the important diagonal terms) by keeping 
A = 1 or another small number, (e.g., A = 4). (Non-diagonal terms may 
need larger ^4-values) 

Rate equation methods based on the density matrix directly work in 
principle with ^N(N + 1) equations (cf. [27] and references), which number 
is for high values of N much larger than that needed in the square root 
method. The frequently used Monte Carlo or "unravelling" method solves 
N x M equations, where M brings in the stochasticity of the environment 
and can be very large , e.g. 100 in [9], 4096 in [43] and 10 7 in [41]). 

(II) A New Form of the Rate Equations. 

The square-root method is formally unlike previously employed methods, 
such as those surveyed in [11], and it can be reasonably expected not to be 
reducible to them. To support this claim, we call attention to the inverse 
"square-root factors" 7 _1 or (7 + ) _1 appearing in equation (13) and other 
equations that follow on. (These factors possess an apparent similarity with 
the "continued fraction method" described in [18] and recently applied in 
[44], but are essentially different.) The presence of these inverse factors finds 
expression also in a novel form of the initial conditions, introduced in section 
2.1. 

(III) Initial time randomization 

As has been shown in section 2.1, the environment induced randomiza- 
tion affects the 7-variables only at the starting time, not during the time 
evolution. This approach was stressed in two early papers by van Hove 
[23, 24] in connection with the derivation of the Pauli master equation. In 
our formulation the initially randomized 7-s lead, in turn, to the ensem- 
ble averaging over various phases in a natural way. Thus the square-root 
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procedure appears to be self consistent. 

In the unravelling method, randomization takes place throughout the 
time development (namely, by the action of random forces on the system). 
It seems to us that the ensemble averaging as done in the square root method 
(namely, by introducing randomness only at the beginning) is closer to the 
meaning originally ascribed to the density matrix. As asserted at various 
places in [23], not only is the repeated stochasticity "superfluous", but "it 
is strictly speaking incompatible with the Schrodinger equation". (We un- 
derstand that van Hove refers here to the form of the Schrodinger equation 
as a first order time- dependent differential equation, so that the system's 
evolution is fully determined by the initial conditions.) 

(IV) Form of the Lindblad operator 

As emphasized at several points in this paper, the choice of the Lindblad 
operator L in the square root method is determined by the requirement 
that the evolution equation for the 7's should lead to rate equations for the 
diagonal terms in the density matrix, that are formally identical to those in 
Lindblad's method [7]. (We assume that this choice is unique.) Thus, we 
do not derive the L-operator from an atomic model for the system and its 
environment. In this sense, the square root method is " phenomenological" . 

However, most density matrix methods are such, too. As an example, we 
call attention to the cases studied in the classic paper of [9]. In some works, 
such as [49] on quantum diffusion, a Lindblad operator was derived from a 
kinematic model, but only through making some approximations. We have 
already noted some dissension, [15]- [22], regarding the employment of the 
Lindblad formalism, especially for initial conditions that are non-separable 
between the system and its environment [19, 22]. 

(V) Pure-to- mixed state transition 

It has been shown in section 2.1 that in the presence of dissipative pro- 
cesses, a pure state becomes a mixed state in a time of order L~ 2 . This is 
also the standard result in other formalisms, obtained theoretically, e.g., in 
[25], and by numerical calculations [27]. 

(VI) Off-diagonal terms in the density matrix. 

For dissipative systems here lies probably the most significant discrep- 
ancy between the square-root and density matrix methods. (In a Hamil- 
tonian system, the two formalisms are equivalent, as already noted.) Our 
defining equations for off-diagonal elements differ from Lindblad's equations, 
and therefore so do (in general) the calculated values of the diagonal terms. 
In a particular case, "level crossing" treated in section 3.1, we have pointed 
out the discrepancy involved in the square-root procedure. Even in this 
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case, the calculated behavior of the diagonal element was very similar to 
that obtained by another method in [27]. 

It seems that the source of the discrepancy is in the different approxima- 
tions made to arrive at a dynamical semi-group from a Hamiltonian system, 
either by truncation of the time domain or by the handling of phase decoher- 
ences [7]. So far we have not found any result that would tend to invalidate 
the present approach, or even to lay bare any shortcomings. In view of the 
questions surrounding the positivity of linear maps with initially entangled 
states [21], the circumstance that our rate equations for non-diagonal terms 
differ from Lindblad's should not invalidate the square root method. 

6 Conclusion 

The "square-root" method (previously used to minimize the action in a gen- 
eral time dependent process [1]) has now been applied to several situations 
where decoherence is expressible by a Lindblad formalism. 

Two irreversible molecular processes (a driven Landau-Zener process and 
the descent to equilibrium across a conical intersection) were formulated 
by the square root formalism. Further illustrations of the method were 
quantum diffusion and memory processes. When compared with results 
given in the literature and obtained using different procedures, the present 
method has led in all cases to good agreement and at the cost of very much 
less effort. The success of the method would seem to justify future uses of 
the formalism, as for Quantum Brownian Motion, already widely studied 
in the literature with density matrix methods, [45] -[53] and regarding the 
question of thermalization [54, 55]. 

At the same time, since the basic equations in the square root method 
and other approaches appear to be different in some respects, there remains 
a theoretical challenge to explain the agreement between the results. An- 
other task is the justification of the basic rate-equation expressions in the 
formalism by starting from microscopic models. 

Appendix 
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A The square-root factored density matrix. 



A.l Matrix formulation of square root factors 

Though the density matrix in the square root factor form has already ap- 
peared before in [2] and [3], we introduce it here following textbook de- 
scriptions of von Neumann's matrix method ([4], [5]). Let \l/ a be a possible 
wave function describing the quantum state of the a'th system in the en- 
semble (a = 1, 2. ..A). It can be expanded in terms of an orthonormal set of 
eigenstates u n as 

N 

*« = 5>n«n (67) 
n=l 

Here the size of the orthonormal set (in principle, infinite) is taken as having 
a finite size N. From the normalization of all the wave functions ^ a we 
obtain: 

N N 

1 = <* a |*a> = (E 7>n\ E T>m) 
n=l m=l 
N N N N N 

= E E «<«nk> = E E 7«<w = E b£l 2 ( 68 ) 

n=l m=l n=l m=l n=l 

As derived in [4] and other texts, the density matrix arises from the ensemble 
average over all systems in the sense that its nm component is 

Pnm = \Y.lX (69) 
a=l 

The 7 Q 's are rectangular matrices of size N x A, distinct for each a (or 
system in the ensemble) and the 7°*'s are conjugate matrices of size A x N. 
Calculating the trace of the matrix p we obtain: 

Tr p= £p nn = f;I^«* = ^El = l (70) 

n=l n=l a=l a=l 

for which we have used equation (68). 

A. 2 Properties of the density matrix 
A. 2.1 The density matrix eigenvalues 

Let us derive some of density matrix well known properties. First p is 
hermitian, since p^ = p, as can be clearly seen from the definition of p given 
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in equation (3) . Thus, all eigenvalues Aj of p are real which is a property 
of all hermitian matrices. 

Second, the sum of all the eigenvalues is unity. This can easily be shown 
as follows. There exists a basis of eigenvectors in which p is diagonal: in 
this basis p will be denoted by pp. The diagonal elements of p£> are just 
the eigenvalues X n . Furthermore, given p in an arbitrary basis there exists 
a unitary matrix U which transforms p to pr, according to the following 
equation: 

PD = Uprf ^p = U^p D U (71) 
By virtue of equation (70) we obtain: 

N 

1 = Tr p = Tr p D U = TV p D = E K (72) 

n=l 

Third, we will show that all X n are positive (or zero). To prove this take 
an arbitrary vector X one can see that: 

N N i A 

xj p x = E E j E Kixx m (73) 

n=l m=l a=l 

In which the definition of p given by equation (3) was utilized. Next define 
R a = J2n=i Xnlni and obtain the result: 

1 A 

x^px = 4 E l^l 2 ^ ( 74 ) 

Further, we write X in the eigenvector basis Vi as X = J2i CiVi which yields: 
X^pX = J2 cMp^CjVj = ]T CIV} ■ ]T XjCjVj (75) 

i 3 i 3 

but since V i ■ Vj = dy we have: 

J2\i\Ci\ 2 = X^pX>0 (76) 

i 

and, since the |Cj|'s are arbitrary, we obtain A« > for every i. From the 
positivity of Aj and equation (72) , we reach the main conclusion of this 
section, that is: 

< Ai < 1 Vi (77) 
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A. 2. 2 Pure and mixed states 

From equation (72) and equation (77) one can reach the following clas- 
sification of density matrices: either there exists a special index s such as 
X s = 1 while for all i / s Aj = or that Aj < 1 for all indices i. The first 
case is denoted as a "pure" state while the second is denoted as a "mixed" 
state. 

For the pure state we have in the diagonal basis p 2 D = pnPD = Pd and 
also in an arbitrary basis 

p 2 = pp= U^p D UU ] PD U = U^p D p D U = p D U = p (78) 

which is a necessary and sufficient condition for a density matrix to describe 
a pure state. One obviously obtains also the result 

Tr p 2 = Tr p = 1 (79) 

For the mixed state we have 

Tr p 2 = Tr p 2 D = ]T X 2 < ]T A, = Tr p = 1 (80) 
In summary, we conclude that 

Tr p 2 < 1 (81) 

in which the equality sign is appropriate only in the pure case. 

As an example for a pure state take an ensemble for which A = 1 and 
pnm = Infm in this case: 

Tr P 2 = S (PP)nn = PnmPmn = 7n7m7m7n 

n n m n m 

= El7n| 2 El7J 2 = l (82) 

n m 

Another obvious case of a pure state is an ensemble with an arbitrary number 
A of wave functions , but in which all the wave functions are equal. Since 
Pnm = \ J2a=i Inlm = Inlm the same argument as above can be applied. 

It remains to show that in case that not all the wave function are equal 
(in a non trivial sense) we obtain Tr p 2 < 1. 



30 



A.2.3 Mixed states 

Let us calculate the trace of the square density matrix: 

Tr P 2 = E E PnmPmn = E E \ E ^ \ E ^ (83) 
n m n m a=l /3=1 

By defining the "state averaged density function" by: 

M^ = E7rf (84) 

we see that: 

A A 

T ^ 2 = 4EEl M ^l 2 (85) 

A a=lf3=l 

Hence we need to show that for all a and (3, M a p < 1 we shall be particular 
interested in the case that the inequality is definite that is M a p < 1. 

Let us look at two arbitrary complex vectors: F, B and a complex scalar 
a = |a|e l<? \ Obviously 

El^ + «^| 2 >0 (86) 

i 

However 

E \Fi + aB t \ 2 = E l^l 2 + M 2 E \ B i\ 2 + a H B ^ F * + a * E B t F i (87) 

i i i i i 

Now denote 

a = El^| 2 >c = El^| 2 (88) 

i i 

which are both real and positive quantities, and 

26 = e ^ E B i F i + E B * F i (89) 

j i 

which is a real quantity. And we arrive at the inequality 

a\a\ 2 + 2b\a\ +c> (90) 

Since, as function of |a|, this is an equation of a parabola which has all its 
values above the \a\ axis (except when the equality holds in this case the 
parabola touches the axis in a single point), it follows that the discriminant 
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of this equation is either negative or zero, the latter in the case that the 
parabola touches the axis in a single point. Thus we obtain 



b 2 < ac (91) 

Now assume that the length of the F and the B vectors is unity. That is 
a = c = 1. Further more denote X = J2i BiF* = \X\e tip . Thus b is equal to 

b= \X\cos((t) + (p) (92) 

The following inequality is obtained 

\X\ 2 < — [ (93) 

By choosing the arbitrary phase 

4> = -(p ± nir (94) 

we obtain the result 

\J2BiF*\ 2 = \X\ 2 < 1 (95) 

i 

this being a special case of the Cauchy-Schwarz inequality. Next let us 
discuss the case in which the equality sign holds in the above equation, that 
is the case in which b 2 = ac and the equation a|a| 2 + 2fe|a| + c = is satisfied 
for a single value of \a\. This can be traced to equation (86) for which we 
have 

J2\Fi + aBi\ 2 = (96) 

i 

but this is only possible for Fj = —aBi. However, since the length of both 
vectors is 1 we arrive at the equation 

Ft = -e**Bi = e^+^Bi (97) 

Thus in order for the equality to hold the vectors F and B must be the 
same up to a "global" phase. We can now show that M a p < 1 and in 
particular M a p < 1, if 7° and 7^ are different in a "non-trivial" way 
(a global phase change does not count as a difference). This is done by 
identifying F = 7", B = 7^. Thus, according to equation (95) , 

|M a/3 | 2 = |£7rfl<i (98) 
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and, in particular, 

|M a/3 | 2 = |E^*l<l (99) 

n 

unless 7® and 7^ are the same up to a global phase. This yields, barring a 
unique case, 

A A 

Tr p 2 = 7iEE i M ^i 2 < 1 ( 10 °) 

71 o;=l/3=l 

Thus a mixed state is really mixed in the sense, that is should contain at 
least two wave functions which are different in a non trivial way. 
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